The effect of short ray trajectories on the scattering statistics of wave chaotic systems 
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In many situations, the statistical properties of wave systems with chaotic classical limits are well- 
described by random matrix theory. However, applications of random matrix theory to scattering 
problems require introduction of system specific information into the statistical model, such as the 
introduction of the average scattering matrix in the Poisson kernel. Here it is shown that the 
average impedance matrix, which also characterizes the system-specific properties, can be expressed 
in terms of classical trajectories that travel between ports and thus can be calculated semiclassically. 
Theoretical results are compared with numerical solutions for a model wave-chaotic system. 



I. INTRODUCTION 

Wave systems appear in many different branches of 
physics, such as quantum mechanics, classical electro- 
magnetism and acoustics. However, solving the wave 
equations in general can be quite difficult, particularly in 
the short wavelength limit for systems which have chaotic 
dynamics in the classical limit i lj. Furthermore, even if 
exact solutions were feasible, there may be uncertainties 
in the locations of boundaries or in parameters specify- 
ing the system. Thus, rather than solving such systems 
exactly, it has often been convenient to create statis- 
tical models which reproduce the generic properties of 
such systems without the need to accurately model the 
details 2]. One successful statistical approach, known 
as random matrix theory, is to replace the exact wave- 
mechanical operators, such as the Hamiltonian or scat- 
tering matrix, with matrices whose elements are assumed 
to be random. Although such formulations cannot pre- 
dict any particular wave system's properties exactly, they 
can predict the distribution of properties in an ensemble 
of related wave-chaotic systems. Random matrix theory 
also predicts the statistical properties of a single wave- 
chaotic system evaluated at different frequencies. The 
random matrix technique applies to a wide range of sys- 
tems and has been well studied both theoretically and 
experimentally. See Refs. [1, 0, H, Q for reviews of the 
theory, history and applications of random matrix theory. 

In this paper, we use random matrix theory to model 
the scattering behavior of an ensemble of wave-chaotic 
systems coupled to the outside world through M dis- 
crete scattering channels. Such scattering systems have 
been studied extensively, with most work focusing on the 
M X M scattering matrix S, either by using a random 
Hamiltonian for the closed system and deriving the re- 
sulting scattering matrix using assumptions for the cou- 
pling between the wave system and the scattering chan- 
nels [3| or by replacingthe scattering matrix with a ran- 
dom matrix directly [1, [^, [l^, [ll|, These two ap- 
proaches are complementary and for some ensembles have 
been explicitly shown to be equivalent Q • 

We consider ensembles of systems whose distribution 
of scattering matrices are well-described by the so-called 



Poisson kernel The Poisson kernel characterizes 

the probability density for observing a particular scatter- 
ing matrix S in terms of the average scattering matrix 
S, which is also called the 'optical scattering matrix'. 
It represents contributions to the scattering matrix from 
elements of the system which are not random. For in- 
stance, if the scattering channels are not perfectly cou- 
pled to the wave system, some fraction of the energy in 
the incident waves will simply bounce off the interface 
between the channel and the scatterer without experi- 
encing the chaotic aspects of the scatterer, thus strongly 
constraining S. This is known as the prompt reflection 
[30] . In addition, rays within the scattering region which 
connect the scattering channels without ergodically ex- 
ploring the chaotic dynamics also affect S [Tl|, [l^ ■ 

Because S is the only parameter in the Poisson kernel, 
methods of finding it for a specific system are of interest. 
Although S can be extracted quite simply from exper- 
imental data, predicting it from first principles is quite 
difficult in general, although it has been done for for some 
specific systems such as quantum graphs In most 

wave systems, however, it depends in a complicated way 
on the interactions between the scattering channels, the 
wave system, and any significant classical trajectories. To 
address this problem, we find it convenient to transform 
from the scattering matrix S to the M x M impedance 
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where Zq is a.n M x M diagonal matrix whose ith diag- 
onal element is determined by the detailed properties of 
the ith scattering channel as described below. We will 
show that the average impedance matrix (to be defined) 
can be expressed directly in terms of classical ray trajec- 
tories. 

Impedance is a meaningful concept for all scattering 
wave systems. In linear electromagnetic systems, it is 
defined via the phasor generalization of Ohm's law as 



V = ZI, 



(2) 



where the M-dimcnsional vector V represents the volt- 
age differences across the attached transmission lines (the 
systems port) and the M-dimensional vector / denotes 
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the currents flowing through the transmission Unes. The 
concept of impedance can be generahzed to cases where 
the cavity is excited through apertures connected to 
waveguides that may support several propagating modes. 
In acoustics, the impedance is the ratio of the sound pres- 
sure to the fluid velocity. A quantum-mechanical quan- 
tity corresponding to impedance is the reaction matrix, 
which is often denoted in the literature as iK J^. 

The diagonal elements of Zq are the characteristic 
impedances of the scattering channels. In electromag- 
netic systems with transmission lines for scattering chan- 
nels, the characteristic impedances are the ratio between 
the voltage difference across the transmission line and the 
current through the transmission line for a monochro- 
matic wave propagating a single direction through the 
transmission line. Other wave systems have analogous 
definitions for Zq determined by the details of the scat- 
tering channels. In what follows, we use terminology ap- 
propriate in the context of an electromagnetic cavity con- 
nected to the outside world via transmission line chan- 
nels. 

With the transformation to impedance, we find that we 
can define an "average" impedance matrix Zavg which is 
related to S via the transformation, 

Za.g^z'J\l + S){l-S)-'z'o^\ (3) 

In contrast to S, we find that Zavg can be evaluated 
directly in the semiclassical limit as a sum over contri- 
butions from the prompt reflection and short classical 
trajectories. Thus, through ([3]) this gives a method for 
approximating S in the semiclassical limit. 

In this paper, we present our approach to calculating 
Zavg- In Secini we provide an overview of our theory for 
lossless systems and describe the most important results 
of our investigations. In Sec. Illli we find expressions for 
the impedance of a specific quasi-2D lossless microwave 
cavity as explicit functions of the boundaries and port 
positions, creating a framework in which we can keep 
some cavity properties fixed and let others change. We 
then apply the semiclassical approximation to our ex- 
act formulations. In Sec. IIII Al we use the semiclassical 
approximation to derive expressions for the impedance 
induced by objects near the port in terms of classical 
short orbits between the ports and the internal scatterers. 
In Sec. IIII Bl we use the semiclassical approximation to 
convert our exact solution with integral operators into a 
finite-dimensional matrix equation with an internal scat- 
tering matrix T. In Sec. IIVI we assume that the matrix 
T is distributed according to the Poisson kernel with the 
average T given by the results of Sec. lIIII and use a result 
by Brouwer Q to find the corresponding distribution for 
the impedance. In Sec. |Vl we demonstrate that in the 
lossless case, our theory agrees with numerical simula- 
tions of our system. In Sec. IVIl we extend our theory to 
lossy cavities and briefly refer to experimental results in 
lossy systems which will be published separately. 



II. OVERVIEW OF LOSSLESS THEORY 

The results presented in this paper are an exten- 
sion of our previously developed random coupling model 
[isl . [TgI . IT^ . For simplicity, we consider only systems and 
frequencies in which the scattering channels have a single 
propagating mode and in which the ports which couple 
the scattering channels to the cavity are separated from 
each other by much more than a wavelength. We found 
previously with these assumptions that by replacing the 
resonant frequencies of our closed chaotic cavity with 
a set of resonant frequencies appropriate to a random 
matrix drawn from the Gaussian orthogonal ensemble 
and modeling the eigenfunctions using the random-plane 
wave hypothesis, the impedance of our lossless wave- 
chaotic systems could be described by [l^l 

Z^tXn + i y/R^io /Rfl , (4) 

where Rn, the radiation resistance, and Xji, the radia- 
tion reactance, are the real and imaginary parts of the 
M X M diagonal radiation impedance matrix Z^, which 
represents the impedance the scattering system would 
have if all the energy which successfully coupled into the 
system was absorbed rather than allowed to couple back 
out. The matrix is an element of the appropriate 
M X M Lorentzian ensemble introduced by Brouwer [3| 
with width 1 and median 0, which in the single-port case 
simplifies to a Lorentzian random variable with width 1 
and median 0. We denote the normalized impedance 
and have previously studied its properties in chaotic sys- 
tems [15, 16, 17, 18]. 

Equation (|4]) is the direct impedance analog of the 
Poisson kernel distribution for S in the case that the only 
contribution to S is due to the 'prompt reflections' caused 
by the impedance mismatch between the scattering chan- 
nels and the wave system [isj . From experimental mea- 
surements and simulation results performed using the 
commercial off-the shelf program High Frequency Struc- 
ture Simulator(HFSS), we know that Eq. ^ describes 
the impedance statistics of our sample systems only if 
the impedances are sampled from a very wide frequency 
range [15]. We find, however, that if we consider sam- 
ple impedances from either narrower frequency ranges or 
from many slightly different chaotic systems, the distri- 
bution of the resulting impedances is still well-described 
by a Lorentzian distribution, but with a different median 
and width than that predicted by Eq. ([4]). 

These deviations are illustrated in Fig. [U which shows 
the median calculated impedance for the quasi-two di- 
mensional cavity illustrated in Fig. [2] (its dimensions are 
given in Sec. |V|. This cavity was the basis for our pre- 
vious numerical and experimental research. It is a sim- 
ulated electromagnetic cavity filled with a uniform loss- 
less dielectric and is coupled to the outside world through 
coaxial cables (the ports) inserted into holes on the top of 
the cavity. For the data in Fig.[l] only port 1 is present. 
On the walls we impose perfect-conductor boundary con- 
ditions. Because our cavity has a uniform height h in the 
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FIG. 1: A comparison between the port radiation reactance as 
measured by HFSS (solid line) and the ensemble median of the 
HFSS simulated impedances (circles). The random coupling 
model predicts that if the ensemble is sufficiently random, 
the ensemble median should equal the radiation reactance. 
The error bars were estimated by assuming that the ensem- 
ble impedance is a Lorentzian random variable (justified by 
statistical examination of the ensemble data) and finding the 
uncertainty in the median, given the numerically found width. 
The differences between these two curves are caused by short 
orbits within the cavity which exist in many realizations of 
the ensemble. 



z-direction which is much smaller than the wavelength of 
the incident microwaves. Maxwell's equations become ef- 
fectively two-dimensional with the electric and magnetic 
fields uniform in the z-direction [l^, [l^. This system 
is an example of a wave billiard, meaning that the rays 
within the cavity follow straight lines except for specular 
reflection at the walls. To produce the simulation data 
shown in Fig. [1] we generated 95 different realizations 
of related systems by adding a small mobile perturber to 
our baseline system and moving it to 95 different, widely- 
spaced locations (see Fig. [2]). We then find that, at each 
frequency, the distribution of impedances is Lorentzian, 
but with a median and width which are different from Xr 
and Rr. (For this example, Xn and Rr are scalars.) We 
also find that as a function of frequency, the fitted me- 
dians and widths oscillate almost symmetrically around 
Xn and Rr. As we will later see, this behavior arises be- 
cause short classical trajectories within the system alter 
the distribution of impedances, analogous to the distor- 
tion of S observed in previous work [ll|, [l^, • 

In this paper we show that corrections to the radiation 
impedance matrix due to the direct orbits redefine Zr in 
an additive way, Zr — > Zavg = .Z^-f [direct orbit terms]. 
More specifically, we find that 

Zavg = Zr + R](^ CRr^ ^ (5) 



where ^ is an M x M dimensionless matrix whose 
(m, n)th element describes the effects of wave propaga- 
tion from port m to port n and is explicitly defined in 
Eq. (|50p . The individual elements Cm,n are each a sum 
over all different possible ray paths going from port m to 
port n, with each path having a phase factor proportional 
to the length of the path and prefactors describing the di- 
rectivity of the classical trajectories relative to the shape 
of the ports, the stability of the ray trajectory and the 
number of reflections from the boundaries. In addition, 
in cases where the geometry of the cavity is varied to 
create an ensemble (for example by moving a perturber 
throughout the cavity) there is a factor that accounts for 
the fraction of realizations in which that particular path 
will contribute (i.e., not be blocked by the perturber). 
We note that as the frequency window within which sam- 
ple impedances are generated gets wider and/or the en- 
semble changes more drastically between realizations, the 
value of C needed to normalize the data goes to and we 
get our original random coupling model back as a limit. 

The substitution of Zavg for Zr can be understood 
to be a generalization of the radiation impedance of the 
ports to include the effects of features of the cavity that 
are distant from the ports but that do not vary from 
one member of the ensemble to another. Consider an en- 
semble of lossless microwave cavities whose generic prop- 
erties (such as volume and circumference) are fixed but 
whose shapes are random and independent except that 
the ports are always placed in the same positions relative 
each other and except that some segments of the wall are 
also fixed. The impedance of each configuration will re- 
flect the interaction between the ports and both the fixed 
and varying segments at the wall. The interactions be- 
tween the ports and fixed wall segments will be the same 
for each member of the ensemble and will contribute to 
the average impedance while the interactions between the 
ports and the varying segments will vary from member to 
member and contribute to statistical deviations from the 
average. As our subsequent analysis will show, the ports 
and fixed wall segments can together be considered to be 
a single super-port which has a radiation impedance of 
Zavg, thus justifying replacing Zr with Zavg hi Eq. 

III. THE IMPEDANCE AS A FUNCTION OF 
CAVITY SHAPE 

For the electromagnetic system described in Sec. [Ill we 
previously derived [l7| the following inhomogeneous wave 
equation for the case where the ports are modeled by ver- 
tical (z-direction), externally imposed, localized current 
densities flowing from the bottom to the top plates, 

M 

(Vi + k^) VT(r) = ikhijY,Mr)Ip, (6) 

p=i 

where is the 2D Laplacian in the {x, y) plane, Vr rep- 
resents the voltage difference between the two plates, 
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FIG. 2: This plot displays a 2-D view of our simulated mi- 
crowave cavity and the perturber positions used to produce 
the ensemble displayed in Fig. [1] The outer walls are fixed in 
all realizations, while every realization has the perturber at 
a different location. The microwaves are fed into the cavity 
through the ports, which are coaxial cables inserted through 
the top of the cavity. The dimension h of the cavity in the 
z-direction (out of the page) are much smaller than the wave- 
lengths used to excite the cavity and therefore results in ef- 
fectively two-dimensional waves in the x-y plane. 



Ip represents the total current injected into the cavity 
through port p, u{f) represents the profile of the current 
injected onto the top plate at port p and has the prop- 
erty J df" u{f') = 1, rj = y/jije is the wave impedance 
of propagation within the medium inside the cavity (i.e. 
rj is the ratio of the electric field to the magnetic field 
in an infinite plane wave), and k = 27r/A is the wave- 
number of the external driving frequency. With perfect- 
conducting boundary conditions, Vr must be zero on the 
cavity boundary. The model of the port considered here 
is appropriate to the case in which the port is smaller 
than a wavelength so that the distribution of current 
(given by u{f')) is fixed, independent of frequency of the 
fields in the cavity Vr {r) . In this way each port is charac- 
terized by a single current Ip, and a corresponding volt- 
age [13, 

Vp=J d'/ Up{/)VT{r'). (7) 

Definition ([7]) was selected in Ref. since it yields P = 
(l/2)Re{V^p*/p} for the power flow into the cavity. The 
cavity impedance then gives the matrix relation between 
the port currents Ip and port voltages Vp. 

Equation ([6]) is the driven Helmholtz equation, and 
although it was derived in the context of quasi-2D elec- 
tromagnetic cavities and a particular port model, it can 
be applied to many different types of systems (such as 
quantum dots or acoustic resonators) simply by relabel- 
ing the constants and tweaking the boundary conditions 



Before considering statistics, we first derive an expres- 
sion for the impedance for individual realizations of the 
cavity. Similar to previous work by Georgeot and Prange 
(20j . we can convert Eq. ^ into an integral equation. 
We do this by introducing the outgoing Green's function 
Go{r, r*; k) which satisfies 

(Vi + fc2) Go(r, r'; fc) = d{f~ r'). (8) 

We then multiply both sides of Eq. © by Go{r, r*; k) and 
integrate both sides over obtaining 

^ dV-'Go(r, r') (W'l + k^) Mr) = (9) 

M 

ikhrj'^ IpJ^d'^f'Go{r,r')up{r'), 
p=i 

where V denotes the two-dimensional domain within the 
cavity. Applying Green's second identity in two dimen- 
sions to the left-hand side of Eq. ([9]) and applying the 
boundary condition on V{r), this becomes 

Mr) = - f dq'Go{r,q')^^^ (10) 
Jav on' 




where q' represents a position on the boundary dV of the 
cavity and the integral over q' integrates over the cavity 
boundary dT>, and where d/dn' denotes a derivative in 
the direction normal to the surface of the cavity at q' . 

In electromagnetic systems, Eq. (|10p has a physical 
interpretation. From Maxwell's equations and the perfect 
conductor boundary conditions, we find that the gradient 
of the voltage, VVr is proportional to the surface current 
in the upper and lower plate of the cavity, with the two 
currents flowing in opposite directions. At the edges of 
the cavity, the surface current flowing in the lower plate 
travels up the outer wall and into the top plate. Thus 
it temporarily travels in the z-direction. We can then 
interpret Eq. pH)) as stating that the electric field inside 
the cavity is simply the sum over the field radiating from 
the changing currents in the ports and the field radiating 
from the changing current in the outer walls. Therefore, 
solving Eq. (fTO| is equivalent to finding the self-consistent 
current induced in the walls, given the currents in the 
ports. We evaluate the normal derivative of Eq. ([T0|) on 
the surface to get our integral equation, 

dVTiq) ^ _ f dGo{q,q')dVT{q') ^^^^ 

On Jqx> dn dn' 

, -I h T f j2^dGQ{r,r') ^ 
+ikh7] > Ip d r Up(r ). 

We henceforth drop the subscripts dT> and T> on the inte- 
gral symbols. As long as the function defining the bound- 
ary of the cavity is well-behaved, Eq. (fTTI) is a Fredholm 
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FIG. 3; A schematic of the different operators defined in 
Eq. p2p . Each operator takes a current at the source po- 
sition and finds the resulting current {K and V+) or vohage 
(Go and V-) induced at the endpoint. 



integral equation of the second type and can be solved 
via the established Fredholm theory [2]| . 

To simply and clari fy o ur results, we follow Prange, 
Fishman and Georgeot |20l . [2^ , and define the operators 



I d-r 



dn 

2^ dGoiq,r') 



dn 



uin, (12) 



V+u{r) 

y_^(g) = VkJ dq'Go{r,q'mq'), 
Gou{r) - J dVGo(r,r')u(r^). 

In these operators, q and q' are real scalars denoting 
points on the cavity boundary dV. They represent dis- 
tance along the boundary of the cavity as measured rel- 
ative to some arbitrary starting position. The vectors f 
and f represent positions within the cavity (i.e. within 
V). Every operator integrates over a primed variable 
and maps it onto the space represented by the unprimed 
variable. 

These operators all have physical meanings in electro- 
magnetism. Go is the two-dimensional, outgoing Green's 
function in empty space; it finds the voltage at some po- 
sition r caused by a delta-function current distribution at 
point 7^. The operator finds the current induced in 
the wall by a delta-function current in the volume. The 
operator K represents the current induced in one part of 
the wall by the current in another part of the wall. The 
operator on the other hand, gives the voltage inside 
the volume which results from the currents in the walls. 
A schematic of the effects of these operators is shown in 
Fig.H 



Using this operator notation and solving Eq. (Ilip . we 
convert Eq. PH)) into 



M 

p=i 



+ 



Go 



Up{r}. 



(13) 

In Fredholm theory, the operator (1 — K) is well- 
defined and can be defined as a ratio of two convergent 
sums 



(14) 



where d„ is an nth order polynomial in the traces of 
If", m < n. For more details on constructing dn, see 
the references [13, [l^ , where it is denoted Dn ■ 

Using the definition of port voltage, Eq. ([7]), we get the 
impedance between ports n and m 



ikhf] 



d^- 



r u„(r) \V_ (l-K) V+ + Go\ u™(f). 

(15) 

The second term in the integral on the right-hand side 
of Eq. represents the impedance the system would 
have if the walls were moved to infinity and outgoing 
boundary conditions were imposed but impedance due 
to direct orbits between the ports were still included. 
Therefore, we define an M x M matrix Zr which has 
the elements 



(16) 



ZR,n,7n = ikhr] / d fM„(r)GoUm(r). 



The diagonal elements of Zfi are equal to the diagonal 
elements of the radiation impedance Zr from Eq. Q [l^ 
and the off-diagonal elements represent contributions to 
the impedance from direct orbits between the ports. Be- 
cause the distance between the ports is large compared to 
a wavelength, we can treat the off-diagonal terms semi- 
classically. The diagonal terms Zf(^n,n depend on near- 
field interactions within the port and thus are sensitive to 
the detailed properties of the port. As in Eq. rather 
than attempting to solve the for the diagonal elements 
of Zji, we treat them as inputs to the theory. This has 
the advantage of freeing us from a detailed port model; 
we expect our model to be accurate even when the port 
behavior is not modeled by Eq. 

In addition, we define the corresponding radiation re- 
sistance and radiation reactance matrices as 



Re 



(17) 
(18) 



where f denotes the conjugate transpose. In the case 
of the Helmholtz equation, or any system in which time- 
reversal symmetry is present, Zr is symmetric and so Rr 
and Xr are simply the real and imaginary components 
of Zfl. 



Equation (jlSp is an exact solution to Eq. ([6]) explic- 
itly in terms of the boundaries. Analytically it is in- 
tractable. In the following two sections, we consider two 
different, but equivalent, approximations to Eq. (jlSp . It 
is by equating these two different formalisms that we de- 
rive our refined theory. 



The short-orbit formulation 



To get useful theoretical results from Eq. p5|) . we make 
the assumption that each port p is located near the po- 
sition ro.p a-nd that is nonzero only within a small 
radius around fo,p. We assume that the ports are phys- 
ically separated from each other and from the walls by 
much more than a wavelength, and that the dimensions 
of the cavity as a whole are much larger than a wave- 
length. With these assumptions, we find that all integrals 
in Eq. (fT5|) (except the diagonal terms in Eq. eval- 
uate Go or its derivatives in the far-field limit. Thus, we 
approximate Go and its derivatives with their asymptotic 
forms, which replaces K with Bogomolny's transfer op- 
erator T [2^, which in the electromagnetic case is given 
by 



T{q,q';k) = -^,/D;Z;e^Sif,f'■,k)-^./4^ 



cos 







(19) 



where 0i{9f) is the angle between the initial(final) wave 
vector and the surface at the position it leaves(hits), 
S{f,r';k) is the classical action along the direct trajec- 
tory from r to r*, and y/Dp p is the stability of the orbit 
from r to 7^ , defined formally as 



Trfc^ 



(20) 



where the derivative with respect to r'±{r'^) denotes the 
gradient with respect to r{f) dotted into a unit vector 
perpendicular to the initial(final) momentum of the clas- 
sical trajectory from r to f* . 

We note that the approximation made in Eq. p9p can 
be used to extend our theory beyond Eq. ([6]) by simply 
changing T to represent the semiclassical approximation 
for other physical situations. For instance, it is possible 
to allow T to have a different action depending on the 
direction of travel, thus violating time-reversal symmetry. 
It is also possible to add loss or gain to the system by 
adding a complex component to the action. 

With these assumptions and approximations, we find 



GoUp{r 



— v/:D^e^^(^>^°-)-'^/4up(-fc.(nr')), 



J21) 

where we have expanded 5 (r, r^) w S{r,fo,p) — ki{r,fo.p)- 
Ar, the quantity ki{r, r*) is the initial wave vector for the 
classical trajectory from rto , and Up{k) is the Fourier 



transform of Up{r) centered on rg.p: 

Up{k)= ( d?Are'''-^^Up{ro.p + Af). (22) 



If we insert Eq. p4|) into Eq. p3|) . we find 



(23) 



Thus the impedance depends only on d„, which de- 
pends on Tr (liT"), and on integrals over the operators 
V—K"-^^V-^-. Evaluating all integrals using stationary 
phase, the prefactors and wave vectors are selected such 
that 1, 23] 



<ffu„ {f)V- K'-^V+u,n (r) 



(24) 



E 

b(/,m,n) 



Un{kf)u*^[k{) 



t/4 



where 6(/, m, n) is an index over all classical trajectories 
that bounce / times, starting at the center of port m and 
ending at the center of port n, S';,(;^m,n) is the action for 
the corresponding classical trajectory and £'b(i,m,n) is the 
stability coefficient defined as in Eq. ((20l) with 5(r, r*) — > 

^b{l,m.n) ■ 

At this point we bring attention to the fact that the 
boundary of the cavity does not need to be connected. 
For example, in the cavity of Fig. [2] there is a circular 
perturber that is moved about, creating an ensemble of 
different cavities. That circle represents a portion of the 
boundary that is not connected to the outer portion of 
the cavity boundary. Equation (|24p in principle includes 
trajectories that pass through the perturber in going from 
one point on the surface to another. (In addition Eq. (IMl) 
includes trajectories that can pass through the convex 
upper boundary of FigJ5)) However, Bogolmony consid- 
ered such trajectories [2J| and found that such unphysical 
orbits terms come in pairs whose semiclassical contribu- 
tions cancel exactly. Thus, the sum over all semiclassical 
bounce terms, which is all that we will consider, will in- 
clude only physical contributions. 

It was found in previous work that the radiation 
resistance for our model ports is given by 



khrj 



2^ 



up{ke)[^ 



(25) 



where 9 is 
(cos(6l),sin(6')), 
gives 



the two-dimensional unit vector 
which when inserted into Eq. 



khr^ j dVM„(f)y_K'-iy+u,„(f) 



(26) 



6(/,m,n) 
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where 



C, 



D. 



b{l,r. 



(27) 



Using similar logic, we can also calculate the off-diagonal 
terms of Eq. (fT6|) . 



where C(o,m,„) and S(^q^^,„ 



(28) 

) are the corresponding pref- 
actor and action for a direct orbit from port m to port n. 
We view the sum over b{l,m,n) as adding successively 
longer length orbits, and we thus refer to Eqs. (I25p - ([28|) 
as the 'short orbit formalism'. 

With this result and similar semiclassical results for 
Tr{K'') [l^l, it is possible in principle to evaluate Eq. pS]) 
semiclassically for any cavity. We will not need to per- 
form this entire calculation explicitly, however. Instead, 
we will use the results from the next section to relate the 
sums over classical trajectories to the elements of random 
matrices. 

We test this short-orbit formalism using the HESS 
program. In the simulator, we construct a fully three- 
dimensional cavity and antenna system similar to the one 
used in previous research; the quasi-2D nature of the cav- 
ity is enforced by choosing the excitation frequency to be 
below the cut-off frequency for modes that vary between 
the top and bottom plates. 

To test our short-orbit theory, we first simulate the 
radiation impedance of a single cylindrically-symmetric 
antenna by placing the antenna inside a circular cavity, 
where the outer circular wall has absorbing boundary 
conditions (to simulate the radiation condition of purely 
outgoing waves) and the port was off-center (this was 
to reduce coherent numerical reflections from the outer 
wall; the numerical absorbing boundary condition is im- 
perfect). As expected, we find a slowly varying function 
of frequency for both the radiation resistance Rr and re- 
actance Xj^. We then change the cavity by introducing 
one perfectly conducting wall into the system, effectively 
producing a cavity in which all waves would either radiate 
away or bounce once off the single wall and then radiate 
away, thus isolating a single term in Eq. (pS)) . Plots of 
such isolated bounce terms are shown in Eigs. [4][6l 

We find empirically that each short orbit term experi- 
ences a frequency-dependent phase shift A(j)p when cou- 
pling through the port p, requiring the introduction of 
a phase factor e*(^'^"+'^'^-) to Eq. 1^. Although this 
phase shift is frequency dependent, in the cylindrically 
symmetric case it is orbit-independent; thus it is possi- 
ble to measure the phase shift using one short orbit and 
then apply it to all others. For our ports, it is most con- 
venient to introduce this phase shift and the cylindrical 
symmetry by simply setting 



>-^6(i,m,n) — c 



b{l.rrL,n) • 



(29) 



This phase shift exists due to the fact that in the HESS 
simulations, we model the ports in detail as a circular 
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FIG. 4: Comparison of the simulated impedance (circles) and 
theoretical impedance (solid line) of a single port with a cir- 
cular perfectly conducting scatterer on one side and radiation 
boundary conditions on all other sides (see the inset in plot 
a). The radius of curvature of the scatterer is 1.02 meters 
and its surface is 7.6 centimeters from the port. The radia- 
tion impedance and phase shift for the port were extracted 
from independent simulation data. Plot (a) shows the re- 
sulting resistances and plot (b) shows the reactances. (Color 
online) 



cross-section coaxial transmission line in which the outer 
conductor contacts the upper plate and the inner conduc- 
tor extends the short way across the cavity and contacts 
the lower plate. The shape and dimensions of port 1 are 
shown in Fig. [T] Port 2, when it is present, has the same 
geometry as port 1, but with an outer radius of 3.0 mm. 
This more detailed port model results in an additional 
phase shift that is not treated in our simple model of 
Eq. ([6]) where we add a fixed current source to the wave 
equation. With this phase correction in Cft(; „i „) , how- 
ever, we find that we can model the impedance very well 
by using Eq. ([^ . 

To show the agreement between theory and simulation, 
we create three different cavity configurations. In the 
first configuration, we have a single port with a single 
conducting, curved wall near it, with absorbing boundary 
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FIG. 5: Comparison of the simulated (circles) and theoretical 
(solid line) ^h,i,2 when the ports are separated by a distance 
of 14.4 centimeters. Note that unlike the diagonal impedance 
matrix elements, the real part of the off-diagonal terms can 
be negative. The radiation impedances and phase shifts were 
extracted from independent simulations of each port. (Color 
online) 



conditions on the remaining surfaces. After simulating 
the isolated ports radiation impedance and finding the 
phase shift, we both predict and simulate the resulting 
impedance in the presence of the curved wall. Fig. 2] 
shows the measured and the predicted impedance from 
5 to 7 GHz, both resistance and reactance, where it can 
be seen that the theoretical and simulated results agree 
well. The deviations between the two are expected; the 
semiclassical approximation is not perfect at frequencies 
this low, and we also have diffraction from the necessary 
truncation of the wall. 

In the second configuration, we tested Eq. ([M)) by in- 
troducing a second antenna into the system and imposing 
radiation boundary conditions on all the outer walls. We 
found the radiation impedance and phase shift of the sec- 
ond antenna using exactly the same methods as for the 
first. We then simulated the mutual impedance between 
the two ports and compared it to Eq. (j28p using the sim- 



FIG. 6: Comparison of the simulated (circles) and theoreti- 
cal (solid line) impedance due to the orbit shown in the inset 
in (a). The ports were 14.4 centimeters apart and 10 cen- 
timeters from the reflecting wall. This impedance is found by 
finding the total impedance of the system with the reflecting 
wall nearby and radiation boundary conditions everywhere 
else and then subtracting the radiation impedance as found 
in Fig.O We denote it ^6,1,2- The radiation impedances and 
phase shifts were extracted for each port from independent 
simulation data. (Color online) 



ulated port parameters, and found that the agreement 
was again excellent, as shown in Fig [51 

The third configuration is the same as the second, ex- 
cept we add a conducting wall next to the two ports, 
creating an orbit which leaves the first port, bounces off 
the wall once, and goes to the second port. To isolate 
the impedance due to this new orbit, we compare the 
changes in impedance (which we denote ^{,^1^2) between 
configurations 2 and 3 rather than the raw impedances. 
Comparing this difference with the semiclassical predic- 
tion from Eq. ([^S)) . we again see excellent agreement, as 
shown in Fig. [B] 

Thus we believe that our short-orbit formulation is ef- 
fective at predicting the impedance of cavities with a few, 
short orbits. In the next section, we discuss an equivalent 
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Side View of Coaxial Transmission Line on Port 1 



h = 7.9 mm 



Inner Radius 0.64 mm 



Outer Radius 2.2 mm 



FIG. 7: Cross-section of the port connecting transmission line 
1 and the microwave cavity. 



formulation which expresses the impedance as elements 
of a finite-dimensional matrix. This equivalence between 
the matrix and semiclassical formulations will allow us 
to create ensembles of cavities which account for short 
orbits within the cavity. 



B. The finite matrix formulation 



The results of Sec. lIII Al are based on an evaluation of 
the continuous, integral operator K defined in Eq. (fT2|) . 
To make connection with random matrix theory, we wish 
to recast the equations in matrix form. The authors of 
Ref. [13 have shown how to do this. Further, in the semi- 
classical limit, the resulting matrix is effectively finite. 

To derive the matrix formulation, we first replace the 
continuous operator K with the semiclassical operator 
T. Then we use the result that semiclassically |22| . 



(30) 



The operator V^^Vlj- was also found in the semiclassical 
approximation to be [2^ 



(31) 



By using Eq. ((30l) to eliminate V- and by adding 
V'|V+ - i {Gq - G]^ /2, which is zero, to the oper- 
ator inside Eq. (|f 5p . we can rewrite the impedance as 



VlV+ = i[Go-Gl 



Zn,7n = \khri J (frun{r) 



+ -I IT' ^ 



(32) 



M Go + G, 



Umir), 

/I + T 
(frun (r ) Vl Y~f y+ (^) ■ 



Fishman, Prange and Georgeot [23] demonstrated that 
in the semiclassical limit, the operator T can be repre- 
sented as an infinite-dimensional matrix whose compo- 
nents are zero except on a finite subspace of dimension 
TV = 2L/A, where in our formulation L is the circumfer- 
ence of the cavity. They demonstrated this by expanding 
all functions on their surface of section (in our formula- 
tion, the cavity boundary) in a Fourier series. In this 
basis, an arbitrary function v{q) is expanded as 



v{q) = ^ a„e 



(33) 



They showed that the operator T, evaluated using sta- 
tionary phase in this basis, is insensitive to Fourier com- 
ponents smaller than a wavelength, resulting in the trun- 
cated subspace. By identical logic, the operator V+ only 
projects onto this semiclassical subspace. Thus, semiclas- 
sically, the function Vl|_Up(r) is non-zero only on this sub- 
space, where it has N discrete components corresponding 
to the Fourier components of the expansion in Eq. (|33p . 

Using Eq. (pij) , we get the dot product between two of 
these vectors, 

y d^fu„(f)y|v+u„(f) = i J cPrUn{r) {Gq - G\^ Um{r) 



khrf 



(34) 



Thus we can rewrite Eq. (??) as 



Z = iXji + V 



1 + T 
I T 



(35) 



where we now treat T as an iV x matrix and where 
D is an X M matrix whose columns Vp are the iV- 
dimensional vectors proportional to the semiclassical 
V-\-Up(f) and which are normalized such that 



(36) 



Because Eqs. ([^5]). (|28l) and (|36p are all evaluated in 
the stationary phase approximation, we can equate the 
matrix elements of any power of T with the semiclassical 
orbit terms 



2vlT 



i-i 



\J RR,nRR 



c, 



b{l,m.n) 



-iTr/4 



(37) 

Equation ([37|) is one of the most important results of 
this paper. By explicitly connecting the semiclassical 
sums to a semiclassical matrix formulation of impedance, 
we can relate the classical trajectories to the more ab- 
stract operator formalism. Thus when we create ensem- 
bles of V and T, we can relate the ensemble averages 
of V and T to the corresponding ensemble averages of 
the classical trajectories within the cavity, which gives 
us a natural method of creating ensembles of T which 
are constrained by short orbits within the system. 
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IV. IMPEDANCE STATISTICS 

As noted in the introduction, it is often difficult to 
solve the wave equation exactly. Even in the semiclassical 
regime, where the problem is in principle tractable using 
classical trajectories, there are difBculties. If the classical 
dynamics is chaotic, the number of classical trajectories 
grows exponentially as does their sensitivity to numeri- 
cal errors. Small mistakes in modeling or small changes 
between similar systems will result in large changes in 
the observed behavior. Thus we follow the long-standing 
tradition of replacing our deterministic expressions with 
statistical models which reproduce the generic behavior 
of the systems being considered. 

Our model for the cavity impedance is given in 
Eq. ()35p . The matrices Xr and v (up to an unmeasur- 
able and thus arbitrary basis, which can thus be absorbed 
into T) are determined by Eqs. p6)) and (fT7|) . which de- 
pend only on the radiation fields from the ports, and 
are thus amenable to direct measurement or non-chaotic 
semiclassical theory. Thus to find our statistical proper- 
ties, we simply seek an appropriate distribution for T. 

The matrix T may be viewed as representing an in- 
ternal scattering matrix. In the case of chaotic dynam- 
ics and in the context of random matrix theory, it is 
most natural to model T as an element of Dyson's circu- 
lar ensemble [s*] (with the time-reversal symmetry deter- 
mined by the symmetry of the underlying system) . If we 
make this substitution, it can be shown that the result- 
ing statistical properties of Z are completely equivalent 
to our previously published random coupling model. In 
this model, the system specific properties of the ports 
(specifically the radiation impedance) are all that is used 
to "normalize" the statistically fluctuating impedance. 
Effectively, the previous model assumes that once wave 
energy enters the cavity it is randomized by the chaotic 
ray trajectories such that no details of the interior of the 
cavity modify the statistics of the impedance. However, 
we have seen in Fig. [T]that there is likely some influence 
of specific ray trajectories within the cavity on the statis- 
tical properties (in the case of Fig. [T]the median) of the 
impedance. We now assume that the effect of these tra- 
jectories can described by the Poisson kernel [l^. That 
is, we assume that the distribution of T is given by the 
Poisson kernel nU 



P{T) 



det 1-TTT 



|^n(/3N+2-/3)/2 



I3N+2-P 



(38) 

where is a normalization constant and T is the aver- 
age value of T over the ensemble, and T is, in principle, 
determined explicitly by the boundaries of the cavities in 
the ensemble. 

The Poisson kernel does not just specify the average 
value of T; it has the general property that [13] 



If we knew T, then we could find the distribution of Z 
directly. Unfortunately, finding T for a specific ensem- 
ble such as that shown in Fig. O is almost as complex 
as finding T and thus has no advantage over numerically 
solving the Hclmholtz equation. By averaging both sides 
of Eq. ([57]) . however, we can find the components of T' 
spanned by the column vectors of v. We find that know- 
ing these average short orbit terms for all I is sufficient 
to get the statistics of Z; because the sum over average 
short orbits is expected to converge, the problem becomes 
tractable for a wide range of ensembles. 

With this assumption for the distribution of T, we can 
find the statistical properties of Z. From a result due to 
Brouwer for matrices distributed according to the Poisson 
kernel we find that we can parameterize T as 



i\HQ+el]W -1 



( Ai^o + el) W + 1 



(40) 



where the scalars A, e and the N x N matrix W are 
ensemble-specific constants which depend only on T, and 
Hq is an TV X iV random matrix distributed according to 
the pdf of the Lorentzian ensemble with median and 
width 1, 



P{Ho) = 



1 



V^det(l + H'2)(/3W+2-/3)/2' 



(41) 



where (3 = 1(2,4) for the Orthogo- 
nal(Unitary,Symplectic) choice of time-reversal behavior. 
Inserting Eq. (gO]) into Eq. we find 

Z = iXR + iX{Wv)^HQWv + ieiWv^Wv. (42) 

We now wish to eliminate A, e and Wv from Eq. (|42p . 
We do this by noting that given the parametrization in 
Eq. (gOI), Brouwer found the value of T to be 



{\ + ie)W''W -1 
(X + ie)W'<W + 1' 



(43) 



Solving Eq. (|43p for M^tj/jA and projecting both sides 
onto the subspace spanned by t>, we get 

(A + le) (Wv)^ ■ (Wv) = Za-og - iXr, (44) 

where we formally define Zavg as 



(45) 



(T') 



(39) 



We denote the Hermitian and anti-Hermitian compo- 
nents of Zavg as Ravg and Xavg, respectively. Matching 
the Hermitian and anti-Hcrmitian components across the 
equality in Eq. ()44|) and noting that Havg must be a non- 
negative matrix because T is subunitary and normal, we 
find that 

X{Wv)^Wv = Ravg (46) 

eiWv)^Wv = Xavg-XR. (47) 
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Inserting these results into Eq. (|42|) gives us 

Z = iXavg + i\/Ravcj^\/Ravg (48) 

where ^ is Hq projected onto the subspace spanned by 
Wv. Brouwer proved that any diagonal submatrix of a 
Lorentzian distributed matrix is also a Lorentzian dis- 
tributed matrix with the same median and width. This 
result combined with the basis invariance of Eq. ((4T|) 
leads to the conclusion that ^ is a Lorentzian random 
matrix with width 1 and median 0, as predicted. The 
basis-invariance of the distribution of also means that 
although Eq. ([^5]) has a family of related solutions for 
Wv, all members of this family are related via a change 
of basis and therefore result in identical statistics for Z. 

Although knowing the form of the distribution for Z 
is useful and can be used fruitfully to fit experimental or 
simulation data, at a single frequency it is only a minor 
improvement over the Poisson kernel in which one can 
also extract S from numerical data [ll|, [l2|- Our last 
step is therefore to predict the value of Zavg using the 
semiclassical approximations developed in Sec. lIII Al We 
do this by noting that because T is subunitary, the mag- 
nitude of all its eigenvalues are less than or equal to one. 
The set of T which have any eigenvalues on the unit cir- 
cle has measure zero. Therefore we can expand Eq. (j45p 
in a convergent series as 

oo 

Za-oa = Zr + 2Y,v^ - T' - v. (49) 
1=1 

Substituting Eqs. ^ and ^ into Eq. (gS]), and re- 
membering Eq. we see that semiclassically Zavg is 
the port impedance plus the sum of the average contribu- 
tions the short orbits make to the impedance. In our case, 
where we have a perturber which moves much more than 
a wavelength between realizations through the entire cav- 
ity, the contributions of orbits reflected off the perturber 
will have an essentially random phase, and we approxi- 
mate their contributions as zero. Because the outer walls 
are fixed, short orbits from the ports to the walls will sys- 
tematically appear in the sum in Eq. (I49p . but must be 
weighted by pb, the fraction of realizations in which they 
do not pass through the perturber. With these results, 
we find the elements of ^ from Eq. ([5]) to be 

6(n,m) 

where the index b(n, m) is over all short orbits which 
go from port m to port n, including direct orbits be- 
tween different ports. Note that when we test this the- 
ory, we use the empirically discovered form of Cft(„ ,„) 
from Eq. (121]). 

We note that Zavg is the impedance the baseline sys- 
tem would have if some fraction of energy were lost every 
time a wave passed through a perturber. Even for very 
large numbers of bounces, this seems to be a general re- 
sult: the impedance needed to normalize the statistics of 



any sufficiently random ensemble will correspond to the 
impedance of a single lossy cavity where loss occurs in 
those features which change between realizations, with 
the degree of loss determined by the degree of change 
in those elements. Thus even with very small pertur- 
bations in which the semiclassical approach is unfeasi- 
ble, the form of Zavg is known, and in analogy to the 
Poisson kernel, we can fit to find the effective radiation 
impedance. More importantly this implies that the fre- 
quency dependence of Zavg matches that of an appro- 
priate lossy cavity. If the perturbations are sufficiently 
uniform within the cavity, the statistics of Zavg evaluated 
over a range of sufficiently separate frequencies would ex- 
hibit the statistics found in our previous work for lossy 
cavities. 

In microwave billiards exhibiting hard chaos with a 
uniform distribution of perturber locations within the 
volume, it is possible to estimate the expected loss pa- 
rameter for Zavg- Because typical ray trajectories within 
the cavity explore the phase space ergodically and be- 
cause the perturber locations are distributed uniformly, 
we expect that, Pb{n,m) ~ exp(— Q;Lh(„^)), where a is 
determined by the perturber size and shape and the 
perturber locations, and L^n,m) is the length of the 
b{n, m)th orbit. Neglecting the phase shifts from the 
traversals of the ports (which are expected to be small 
due to the small size of the ports), we find that with this 
expression for Pb(n.m), 

(n,ni^ ^ Cf,(„,™) exp(i(fc-|-ia)Lf,(«,m) -«7r/4), (51) 

6(n,m) 

where we have used that the classical action in microwave 
billiards is given by Sb{n,m) = kLi,i^n,m)- We found in 
previous work [l6l| that the transformation k ^ k + ia 
corresponds to adding loss, with the loss parameter given 
by Q = k/{2a) and the line- width to level spacing ratio 
a = kAa/n, where A is the area of the microwave cavity. 
Thus, by analogy, introducing pi, into Eq. ([50]) is roughly 
equivalent to Zavg being the impedance of a lossy cavity 
with uniform loss. Using a Monte Carlo simulation of 
long ray trajectories in our bowtie billiard, we find that 
for the perturber positions in Fig. [21 ~ .25m~^. 

This result for the average impedance is related to work 
done by Brouwer and Beenakker [2^ following Biittiker 
[26| . In their work, they found that a lossy quantum 
dot could be modeled as a lossless quantum dot cou- 
pled both to the physical scattering channels and to a 
large number of weakly coupled parasitic channels. Be- 
cause Zavg represents the impedance of a system where 
a small amount of energy is lost when it passes through 
a perturber position, the perturber positions function as 
effective parasitic channels, and we expect the resulting 
average impedance and average scattering matrix to be 
well-fit by their theory. 
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V. NUMERICAL TESTS OF THE 
THEORETICAL PREDICTIONS 

In this section, we show the results of several tests 
of Eqs. (O and ([50)1 for both one and two-port con- 
figurations. We obtained the data for these tests from 
simulations using HFSS. The cavity configuration we 
used is shown in Fig. [21 Port 1 is centered at loca- 
tion {x — 18.03cm, y — 15.48cm) and port 2 (when it 
is present) is at (x = 36.7cm,?/ = 15.48cm). The lower- 
left corner of the cavity is at (x = 0.0cm, y = 0.0cm). 
Both ports have essentially the geometry as shown in 
Fig. 1 of Reference [l6[, but with different dimensions. 
Both ports have an inner radius of .635mm, but port 1 
has an outer radius of 2.29mm while port 2 has an outer 
radius of 3.05mm. The lower and left straight sides of 
the cavity have lengths Li = 43.18cm and L2 — 21.59cm 
respectively and the upper and right sides have radius 
of curvature Ri = 103cm and R2 = 63.9cm respectively. 
By moving a perfectly conducting circular perturber with 
a diameter of 2.54cm to 95 different locations (shown as 
circles in Fig. ^ within the cavity, we construct our en- 
semble. The resulting impedances were simulated at 201 
uniformly spaced frequencies from 5 to 7 GHz, inclusive. 
For this frequency range, we get that the effective loss 
parameter Q due to the ensemble averaging changes lin- 
early from about 210 at 5 GHz to 300 at 7 GHz. The 
effective line-width to level spacing ratio a also increases 
linearly from about 0.95 to 1.35. 

A. Single-port Tests 

To test our predictions in the single-port case, we first 
confirm that at each frequency the impedances have a 
Lorentzian distribution, fitting to find the median and 
width. From the fit, we see that Eq. ((50)) does not con- 
verge quickly enough to be practical at a single frequency. 
However, we find that if we use frequency averaging and 
short orbits together, we regain universal statistics over 
a much narrower frequency band than was required with 
out previous theory [l^ (i.e. if Zr is used as in Eq. Q). 
In addition, we confirm that as a function of frequency, 
the fitted Zavg has the characteristic behavior of the 
impedance of a lossy cavity. 

To test that the impedance is Lorentzian distributed 
at each frequency, we numerically find the three quartiles 
Q(i,2,3)(/) of the 95 sample impedances (denoted Zi{f)) 
and thus find the sample median Zmedif) = Q2{f) and 
the sample width Z^id{f) = Qsif) - <9i(/)- Assuming 
that Zmedif) and Zwid{f) are approximately the correct 
median and width, we can then find the phase of the 
normalized scattering matrix coefficient 

</>.(/) = 2tan-i(^^i:^^-^^). (52) 

For large N and with a Lorenzian distribution for Z{f), 
4>i{f) should be uniformly distributed between — tt and 



TT. We bin the numerical (ja into 10 equal-sized bins 
and find the resulting deviation from a uniform dis- 
tribution. We then find the anticipated distribution of 
by performing a Monte-Carlo statistical analysis on 
402000 realizations of 95 Lorentzian random variables 
each, transforming them precisely as done in Eq. (j52p 
and finding the resulting . Comparing the values of 
for the simulation data and the Monte-Carlo distribu- 
tion, we find that at the 95% confidence level, we can ac- 
cept the Lorentzian hypothesis for 93% of our frequency 
impedance samples, while at the 99% confidence level, 
we can accept all of our impedance samples. These re- 
sults are consistent with the data being distributed with 
a Lorentzian distribution at each frequency. 

To test the semiclassical theory in the single-port case, 
rather than using the sample median and width we at- 
tempt to predict the median and width using Eqs. ([S]) 
and (j50p . In practice, we must eventually truncate the 
sum over semiclassical trajectories. We therefore define 
the truncated average impedance 

Zt.i,, = ZR + 2Y,A-T' -vi, (53) 

which is equivalent to the sum over all classical trajecto- 
ries which bounce up to TVfc times. Thus the normalized 
scattering phase we use to test our semiclassical theory 
is 

0i,Ar, (/ ) = 2 tan (-t ^ I, (54) 

where Rt.N^ and Xt.N^ are the real and imaginary parts 

of Zt^N^. 

In Fig. [H we compare Zt^N^ to Z^iid for Ni, = 2, 5, 6. 
We see that as Nb increases, the two terms becomes in- 
creasingly similar. Zt,2 follows the frequency average of 
Zwid- Also, despite the fact that Z^id changes rapidly 
in frequency, consistent with the assumption that it is 
the real part of a lossy impedance, Zt.5 and Zt^ have 
begun to fit even the large spikes in Z^iid- This strongly 
supports the validity of Eq. (|50p and implies that with 
a sufficiently large number of bounces or a more random 
ensemble (such as including more perturbers or having 
a larger perturber), it may be possible to predict Zavg 
at a single frequency. We have not yet confirmed this 
possibility. 

In previous work, we found that normalizing the 
impedance with zero bounces (i.e. the radiation 
impedance Zr) was sufficient to get universal statistics 
if we sampled the impedances over a sufficiently wide fre- 
quency range. If the frequency window was too narrow, 
however, we found systematic deviations from universal 
statistics. Even though we do not have enough terms in 
Ztfi to experimentally find Z^vg semiclassically at a sin- 
gle frequency, we find that by combining short orbits and 
frequency averaging, we can get universal statistics over 
much narrower frequency ranges than previously. 
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FIG. 8: A comparison between the fitted widtli of tlie 
impedance distribution Z^id at each measurement frequency 
and the semiclassical predictions for the widths due to short 
orbits which bounce up to 2, 5 or 6 times. By 6 bounces, 
the semiclassical prediction has begun to fit the gross fea- 
tures of the fitted widths but is far short of the number of 
orbits needed to fit the sharp spikes. In addition, Gibbs phe- 
nomenon has become a problem in the sixth bounce around 
5.75 and 6.7 GHz with the semiclassical prediction dipping 
too close to zero. (Color online) 



FIG. 9: The average deviation of the (/>i,jVb from the uni- 
form distribution for different frequency windows and differ- 
ent choices of N},. The were calculated for all possible fre- 
quency windows with a given width and then averaged over all 
different window realizations. The random control was gen- 
erated by a Monte-Carlo simulation of the x for uniformly 
distributed phases. Note that for the sixth bounce we exclude 
frequency windows in which the calculated Ravg falls below 
0.1. (Color online) 



To measure the deviation of the distribution of the 
measured (/!'i,Ar„ from uniform, we introduce the statis- 
tic. For a frequency window of width (5/ and centered at 
/o, the statistics is calculated by binning the <i>i^Ni, 
from every realization and from every other frequency in 
the frequency window into ten equally sized bins from — tt 
to TT. (We take every other frequency window because the 
impedance values of adjacent frequencies are found to be 
strongly correlated.) The x^ statistic for this window is 
then given by 

r— 1 ^ ' 

where N^. is the number of in the rth bin, and (A^, ) 
is the expected value of given a uniform distribution. 
The x^ statistic is chosen because it has approximately 
the same distribution independent of A/. 

In Fig. m we display the average x^ statistic from our 
sample for multiple values of l^\y and different choices of 
window width. The averaging is performed over all fre- 
quency windows of the same width, including windows 
whose frequencies overlap. We see that for small window 
widths, increasing iVf, systematically decreases the error. 
In addition, we also see that, up to a point, increasing 
the window size also decreases the error, but, once the 
error has leveled off, it decreases no further. This is con- 
sistent with the frequency averaging effectively removing 



the longer orbits, making the improved statistics from 
the larger N\, irrelevant. In addition, for comparison we 
include the x^ statistics for a set of truly independent 
random phases and see that for the largest and widest 
window widths, our results are statistically indistinguish- 
able from true randomness. 

At this point we note one caveat to our use of Eq. ([M]) 
relevant to Fig. [5] In Fig. [51 we see that Zt^% drops almost 
to zero near 5.7GIIz and 6.7GHz. In fact, if we continue 
to add bounces, Rt,Nb can actually become negative at 
some frequencies, which would represent gain in a lossy 
system and is unphysical. This unphysical behavior oc- 
curs only because the sum in Eq. ([5(7)1 has been trun- 
cated, but it badly distorts the calculation in Eq. ([M|) 
within the affected frequency ranges due to the abnor- 
mally small denominator. This occurs because our sum 
over orbits is effectively an attempt to expand a function 
with poles near the real axis in a Fourier series. Due to 
the rapidly changing features in Ravg,i,i, we get a form of 
Gibbs phenomenon, in which a Fourier series attempting 
to fit a discontinuous function systematically overshoots 
the fitted function. For the purposes of producing Fig. [31 
we simply ignored frequency windows which contained 
frequencies such that Rt,e{f) < 0.1i?fl(/). It may be 
possible to avoid this problem by using a smarter method 
to expand Eq. ([^ . 
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FIG. 10: A comparison between the numerically found me- 
dian (solid lines) and width (dashed lines) of the distribution 
of impedances at each frequency for the single-port simulated 
impedance (the black lines) and the (1,1) element of the two- 
port simulated impedance (the gray lines). Note that for clar- 
ity we have subtracted the radiation impedance of the single 
port from the medians of both sets of data. 



B. Two-port tests 

To test the results from the two-port configuration, we 
first test the statistical properties of the diagonal ele- 
ments of the two-port impedance, considered separately. 
From Eqs. (|48l) and (|36p . and remembering that the dis- 
tribution of ^ is basis-independent, we find that the diag- 
onal elements of the multi-port impedance, considered in- 
dependently, should be Lorentian random variables with 
width Ravg,n,n and median Xavg,n,n- Thus the diagonal 
elements of Z are susceptible to the same analysis used 
in the single-port case. When we perform this statistical 
analysis on both diagonal elements of Z considered in- 
dependently, we get results essentially identical to those 
shown for the single-port case. 

Because port 1 is in the same location for both the 
single-port and two-port simulations, our theory also pre- 
dicts that the median and width of the distribution of 
Zi i at a single frequency should be almost identical to 
the median and width of the single-port impedance distri- 
bution at the same frequency. The data is in agreement 
with this prediction as shown in Fig. [10] (note that, for 
the sake of clarity, we have subtracted the radiation re- 
actance of the single port from the plotted medians of 
both sets of data). 

The statistics of the two-port normalized impedance 
are more than the independent statistics of the diagonal 
elements; the elements of the 2x2 matrix Z are strongly 
correlated. All elements of Z go to infinity, for instance, 
when the frequency goes through a cavity resonance. One 
common way of expressing this is via the correlation be- 
tween the eigenvalues of ^. The eigenvalues of ^ have the 



form cot(6'„/2), where the distribution of the 0„ is given 

by 



P{{9n}) oc n - e 



(56) 



For the two-port impedance problem, this distribution 
simplifies to 



1 



sm 



71 - U2 



(57) 



Thus to test our theory, we must fit or calculate Zavg, 
find the values of ^ for our sample data, diagonalize, and 
find the distribution of the differences between the result- 
ing phases. We again use to determine the goodness- 
of-fit with the definition from Eq. ([55]) . but with {N,.) 
determined by integrating Eq. ([57)1 . As we did for the 
single-port case, we both fit to find the numerical Zavg 
and use the semiclassical sum. 

Numerically fitting Zavg is more difficult for the two- 
port case than in the one-port case. We can find the 
diagonal elements of the fitted Zavg simply by fitting the 
diagonal elements of Z to Lorenzians exactly as in the 
single-port case. Fitting the off-diagonal elements nu- 
merically is more complex because both the shape and 
width of the distribution of the off-diagonal elements of 
Z depend in a non- linear way on all the elements of Ravg ■ 
Rather than attempting this more complex fit, we con- 
sider the rotated impedance matrix 



ozo' 



lOXavnO' 



(58) 



j (o v^o^) (oio^ ) (o /r::;;o^) , 



where O is a constant orthogonal matrix. Because the 
statistics of ^ are independent of basis, the diagonal el- 
ements of the rotated matrix OZO^ will be Lorentzian 
distributed random variables with widths and means give 
by the diagonal elements of OXavgO'^ and ORavgO^ . 
Thus if we make the simple choice 



O 




(59) 



we find that the widths of the distribution of the diagonal 
elements of the rotated impedance are 



■71,1 



R. 



■avg.1,1 



■avg,2,2 



{ — y^ Ravg,l,2j 

(60) 

where n — 1,2, and with corresponding logic for Xavg- 
Thus to find Zavg,i.2, we fit the diagonal elements of 
OZO^ to Lorenzians and take the half difference be- 
tween the fitted medians and widths for the different di- 
agonal terms. This algorithm generalizes to larger num- 
bers of ports by choosing O to rotate between the ap- 
propriate port indices. In the process, we also confirm 
that the diagonal elements of the rotated impedance are 
in fact Lorentzian distributed, providing further support 
for Eq. 
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FIG. 11: A comparison between the universality of the phases 
of ^ as extracted using different methods for normalizing the 
two-port impedance data. All data shown is the statis- 
tics for fitting the different between phases to Eq. (|57|) . The 
gray lines represent error of the data normalized using the 
semiclassical impedance sum. The solid line represents the 
phase difference statistics, but at each frequency using the 
fitted value for Zavg, where the fitting parameters are found 
as described in the text. (Color online) 



In Fig.[TTl we show the results of this analysis. Rather 
than considering frequency windows, we simply consider 
the statistics of 16*1 — 02 1 at each individual frequency us- 
ing different values of Zavg ■ We find that the numerically 
fitted Zavg normalizes the data well, but not perfectly, 
with the range of values falling well within accept- 
able bounds as determined by the theoretical distribution 
of x^, but systematically larger than would be expected 
from true randomness. For the truncated sums, how- 
ever, the story is rather different. We see that at many 
frequencies, no corrections are needed at all to get good 
statistics. Some frequencies, however, have large devia- 
tions from Eq. ([57)) . That these deviations are caused 
by short orbits is demonstrated by the fact that, as we 
add 1 and 2 bounces to the semiclassical Zavg, the devi- 
ations initially become smaller, almost reaching the level 
of noise. Unfortunately, adding longer orbits does not 
improve the situation; they leave the statistics either un- 
changed or markedly worse. The reasons for this are 
unclear but are likely due to a combination of the Gibbs 
phenomenon observed in Sec. IV Al and the tendency of 
some impedance matrices to be poorly conditioned due 
to systematically large values of cot(6'„/2). 

Thus we have confirmed many of the predictions of 
the extended random coupling model, including the ap- 
proximate independence of the statistics of the diagonal 
elements of Z , the invariance of the distribution under 
rotation, and the level spacing statistics for pairs of eigen- 
values of ^. 



VI. ADDING LOSS 

In practice, no real cavity will be truly lossless. For 
non-zero frequencies, even cavities with superconduct- 
ing boundaries have loss due to interactions between mi- 
crowave photons and quasiparticles. In quantum me- 
chanical systems, there will always be dephasing, which 
is functionally equivalent to loss [1^. In this section we 
therefore address the effects loss has on our theory. 

From Maxwell's equations, we derived in previous work 
[l^ that going from lossless to uniformly lossy is per- 
formed by the transformation k ^ k + ia, where a — 
kj (2(3) and Q ^ 1 is the loss parameter of the closed cav- 
ity. Performing this analytical continuation takes some 
care. The function that must be explicitly continued ana- 
lytically to obtain universal statistics of the sort we found 
previously is the normalized impedance ^, given by 



^^avg 



[Z - iXa 



\R, 



-1/2 



(61) 



Because Ravg and Xavg appear independently in 
Eq. (|61|) . we must analytically continue each indepen- 
dently. Because taking the real and imaginary parts of 
non-constant functions is not an analytic operation, for 
lossy systems Ravg and Xavg are not real but rather the 
analytic continuation of the real and imaginary part of 
the lossless Zavg on the real axis. These analytic contin- 
uations are unique. 

To find this analytic continuation for the microwave 
billiards used in our experiments, it is necessary to ex- 
plicitly find the real and imaginary parts of Zavg- For 
our microwave billiard system, we find that 



R, 

X, 



■avg 



avg 



— Rr ■ 

= Xr 



(62) 
(63) 



where the M x M matrices p and x have the elements 



Pn,; 



Xn,: 



Pb{n,rn)y Db(n,m) 

X cos (fc (L 

Pb{n,m) ^b{n,m) 

b{n,ni) 

X sin {k {Lp,n + -^p,m + ^6{n,m)) " 



(64) 



^/4) 



(65) 



7r/4), 



where we have made the empirically observed substi- 
tution A0„ = kLp^n, Lp^n is observed to be a port- 
dependent constant, and L^n,m) is the length of the tra- 
jectory b(n,m), where we have used that for billiards. 

Because Zn and Di,(n,m) change slowly in frequency 
compared to the level spacing, they are approximately 
independent of a and thus equal to their lossless coun- 
terparts. Thus the analytic continuations of Ravg and 
Xavg consist of keeping the forms of Eqs. (|62p - (l5^ un- 
changed but allowing k to become complex. 
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After performing this continuation, ^ will no longer be 
real. However the distributions of the real and imagi- 
nary parts oii£ have been found as a function of the loss 
parameter, [271 [28j and, for low-loss systems, we expect 
the distribution to be approximately universal. The only 
difficulty with this analytic continuation is that the sum 
in Eq. (|49|) will not necessarily converge if the loss pa- 
rameter is too high. In such a case, it is necessary to 
perform the analytical continuation on the form of Zavg 
given in Eq. ((45|) . At this time, we have not attempted 
this, but we anticipate that it will require evaluating the 
denominator in Eq. (j45p via a method other than short 
orbits, a poosible subject of further research. 

We have experimentally tested this theory in the one- 
port case for a microwave quarter-bowtie billiard de- 
scribed in previous work for frequencies from 6 to 18 GHz 
[29| . The primary difhculty in applying Eqs. (O and (l50l) 
over this frequency range is the fact that the loss param- 
eter is not constant as a function of frequency. However, 
our theory predicts that the phase of the normalized scat- 
tering parameter s will be uniformly distributed indepen- 
dent of the loss parameter [l3| , where s is given by 



Thus by fitting the distribution of the phase of the nor- 
malized scattering parameter to the uniform distribution, 
we can again find the statistic for various frequency 
windows and different choices for the number of bounces 
before truncation Nh. We have performed this experi- 
ment with a lossy cavity and found results qualitatively 
similar to Fig.[9][2i|. 



VII. CONCLUSIONS 



In this paper, we have shown that the random cou- 
pling model, Ref. [13] and Eq. (g]), can be extended to 
take into account system specific short orbits that af- 
fect the statistical features of the system. From the 
numerically and experimentally observed deviations of 
our results from universality, we anticipated that inter- 
actions between the walls and the port that were not 
sufficiently changed from realization to realization would 
result in corrections to our model. We then derived a 
model that could predict such corrections. Numerically 
and experimentally, we found that the improved model 
resulted in statistically significant improvement when fit- 
ting to the random coupling model, effectively reducing 
the deviations to the level of noise. In addition, we devel- 
oped utilizations of several mathematical tools, including 
Prange's semiclassical version of Bogomolny's T opera- 
tor, that could be fruitful in further study of chaotic cav- 
ities and wave-chaotic systems with known dynamics in 
general. 
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